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Recent experiments have confirmed the importance of nuclear quantum effects even in large 
biomolecules at physiological temperature. Here we describe how the path integral formalism can 
be used to describe rigorously the nuclear quantum effects on equilibrium and kinetic properties of 
molecules. Specifically, we explain how path integrals can be employed to evaluate the equilibrium 
(EIE) and kinetic (KIE) isotope effects, and the temperature dependence of the rate constant. The 
methodology is applied to the [1,5] sigmatropic hydrogen shift in pentadiene. Both the KIE and the 
temperature dependence of the rate constant confirm the importance of tunneling and other nuclear 
quantum effects as well as of the anharmonicity of the potential energy surface. Moreover, previous 
results on the KIE were improved by using a combination of a high level electronic structure cal- 
culation within the harmonic approximation with a path integral anharmonicity correction using a 
lower level method. 



I. INTRODUCTION 

After solving the electronic structure problem, most 
molecular modeling simulations treat nuclei as classical 
particles. While this often is an appropriate point of 
view, it fails completely in some cases, and almost al- 
ways when hydrogen is involved in bond breaking or for- 
mation. Strong nuclear quantum effects have been ob- 
served in many reactions and recently even in several 
enzymatic reactions [ll-O- In such situations, one should 
treat some or all nuclei quantum mechanically. In partic- 
ular, the nuclear quantum effects should be considered in 
the simulation even though they may be similar for the 
corresponding reactions in enzyme and in solution, and 
therefore may not contribute to the catalytic effect of the 
enzyme 0, 

Solving the time-dependent Schrodinger equation ex- 
actly is of course possible only for a handful of atoms 
Q. Luckily, in most chemical reactions, some quan- 
tum effects are at least partially washed out by finite 
temperature. Whereas the "real-time" quantum dynam- 
ics is extremely difficult, quantum thermodynamics, or 
"imaginary-time" quantum dynamics, can be computed 
accurately for fairly large systems using the Feynman 
path integrals 0, llj ■ 

In general physics, path integrals are mostly known as 
an elegant tool to formulate analytical theories in particle 
physics and quantum field theory. In chemical physics, 
the somewhat abstract path integrals have evolved into 
an extremely practical numerical tool [l^. In fact, 
they are probably the most successful tool in solving 
quantum statistical problems in large systems without 
symmetry. 
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In this paper, we describe how path integrals can be 
used to compute the equilibrium isotope effects, kinetic 
isotope effects, and the temperature dependence of the 
reaction rate constant. 



II. METHODOLOGY 

A. Path integral methods 

The central quantity in quantum statistical mechanics 
is the partition function. 



(1) 



where /3 = (/csT)""^ is the inverse temperature and £„ is 
the energy of the system in the eigenstate n. If the parti- 
tion function is known analytically, any thermodynamic 
quantity can be found. E.g., the thermal energy can be 
computed as 



E = - 



d(3 ' 
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The beauty of path integrals lies in that they allow 
computing the partition function without finding the 
eigenstates of the Hamiltonian. Let us therefore con- 
sider a molecular system consisting of N atoms with 
masses m^. Starting from the exact expression Q(/3) = 
Tr(e~'^^), one obtains the path integral (PI) representa- 
tion of Q as Q = limp_i.oo Qp, with [lOj 



Qp(/3) = C / dr(i) 
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i-(s) = (r^^*) ^ r^'*) ^ _ , ^ ) the set of Cartesian coordinates 
associated with the sth time shce. Finally, $ ({r'^'*^}) is 
the effective potential given by 
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with r*^°) = r'^^'' and {r'^''^} representing 
{r(i),r(2),...,r(-P)}. 

From expressions (jSj and (HJ it is obvious that for 
P = 1 one obtains the classical partition function and 
therefore classical thermodynamics. The quantum ther- 
modynamics is obtained in the limit P — >■ cx), but in 
practice it often suffices to take a finite value of P to 
obtain accurate quantum results. 

To compute a thermal average A{P) of a physical ob- 
servable A such as the heat capacity or a rate con- 
stant, one starts from the exact quantum expression 
A{I3) = Tr[Aexp(-/3i?)]/(5(;3) or, if possible, from an 
exact expression for A[j3) in terms of the partition func- 
tion, such as Eq. ^ for energy. Using the PI expression 
for Q, one ends up with a PI expression for the ther- 
mal average 
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where {■) ^ denotes a weighted average over PI integral 
configurations. The weight is given by p = exp(— /3<I>) 
and the quantity yl({r(*'}) is called a PI "estimator" 
for A{I3). 

The PI average ([5]) can be evaluated efficiently using 
a PI molecular dynamics (PIMD) or PI Monte Carlo 
(PIMC) techniques In the calculations presented 

below, a PIMD implementation in the molecular dynam- 
ics package AMBER 10 [llj was used for the equilibrium 
isotope effects while an in-house PIMC code was used 
for the kinetic isotope effects and the temperature de- 
pendence of the rate constant. 

Up to this point, everything was straightforward. The 
interesting twist comes because the estimator A ({r^**^}) 
is not a unique function of {r^**)}. The art of Pis lies in 
finding the optimal estimator in the sense of having the 
smallest statistical error for a given number M of sam- 
ples. As the statistical error of a simulation is propor- 
tional to M~^/^, an estimator with a smaller statistical 
error can lead to a much more efficient simulation. 



B. Equilibrium isotope effects 

The first application that we shall consider is the cal- 
culation of the equilibrium isotope effect (EIE), defined 
as the effect of isotopic substitution on the equilibrium 



constant. More precisely, 

EIE : 



(6) 



where K is the equilibrium constant and I (h) denotes 
the lighter (heavier) isotope. The equilibrium constant 
can be expressed in terms of partition functions as K = 
Qp/Qr where subscripts p and r denote the product and 
reactant, respectively. Hence the EIE can be computed 



as 



EIE : 



exp 



exp 
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where F = — /3~^log(5 is the free energy and A is a 
parameter providing a smooth transition between the 
lighter and heavier isotopologs [l3|. For instance, it can 
be accomplished by a linear interpolation of masses of all 
atoms in a molecule according to the equation 



■m^{X) = (1 - A) rrih^i + Xmi^, 



(8) 



Unlike the partition function or the equilibrium constant 
K, the derivative dF/dX of the free energy can be com- 
puted by a PIMD or PIMC simulation directly, and so the 
equilibrium constant can be computed by the so-called 
thermodynamic integration (TI) expressed in Eq. ([7|). 
The estimator for dF/ dX can be found in Ref. [T2I . 



C. Quantum instanton approximation 

Theoretically more challenging quantities than the 
equilibrium constant or EIE are the rate constant or the 
kinetic isotope effect (KIE). The reason is that these ki- 
netic quantities combine quantum thermodynamics with 
the real-time quantum dynamics. However, since the 
most interesting chemical reactions occur at a thermo- 
dynamically relatively high temperature, there exist suit- 
able approximations. 

A recent and accurate approximation which takes into 
account most quantum effects is the so-called quantum 
instanton (QI) approximation of Miller and coworkers 
p^ . In this approximation, the rate constant is ex- 
pressed as 



fcQi(/3) 



A h Cff(o) 

2 AH Qr ' 



(9) 



where Cff(t) is the flux-flux correlation function, Qr the 
reactant partition function, and AH a certain energy 
variance near the transition state. 



D. Kinetic isotope effects 

The kinetic isotope effect is defined as the effect of 
isotopic substitution on the rate constant, 

k. 



KIE 



kh' 
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This KIE is widely used in chemical kinetics to detect nu- 
clear tunneling and other quantum effects, as well as to 
distinguish between possible reaction mechanisms. Un- 
like the rate constant itself, the KIE depends very little 
on the classical energy barrier height, and so the KIE 
can separate various effects from the simple exponential 
dependence on the barrier that overwhelms the rate con- 
stant. Using the QI expression the KIE can be ex- 
pressed as [IJ] 



KIE = 



Cff.i(O) 

Qr,h Cdd,/(0) Cdd,l(0) 



Odd,h(0) 



(11) 



where the delta-delta correlation function Cdd(i) was in- 
troduced. This correlation function at time < = is a gen- 
eralization of the partition function, constrained to two 
dividing surfaces. Introduction of Cdd into expression 
(llip simplifies the calculation: On one hand, similarly to 
AH the flux factor Cff(0)/Cdd(0) can be computed di- 
rectly in a sing le PIMC or PIMD simulation 15 . On the 
other hand, the ratio Cdd,;(0)/Cdd,h(0) can be computed 
by a TI analogous to Eq. ([7]). The estimators can be 
found in Refs. M, M, [H 

It should be noted that several other PI approaches 
to compute the KIE exist, which are however not based 
on the QI. These include, e.g., approaches based on other 
quantum transition state theories fi g. 19,1 or on the quan- 
tized classical path method l20l l2l|. 



between the reactant and transition state energies and 
the logarithmic derivatives of Qr and Cdd, one can com- 
pute the ratios of Qs and CddS as 
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The estimators for the reactant energy Er are well-known 
for PIMD simulations |23-l26l|. Estimators for the tran- 
sition state energy E^ were developed in Ref. [13 for sev- 
eral types of constraints. These estimators are derived in 
a similar manner, based on the rescaling of coordinates 
and a finite difference evaluation of the derivative with 
respect to /3 that was used for the first time by Predescu 
and coworkers for reactant energies and heat capacities 
p7l [. The statistical error of these estimators showed a 
surprising behavior in comparison with the estimators for 
Er- In particular, the centroid virial estimators was not 
always the optimal estimator [13l . Below we use estima- 
tors from Ref. '22' for both Er and E^ because they are 
suitable for a PIMC simulation. 



III. COMPUTATIONAL DETAILS 



A. Equilibrium isotope effects 



E. Temperature dependence of the rate constant 

The last application that we shall consider is the di- 
rect evaluation of the temperature dependence of the rate 
constant, i.e., of the ratio fc(T)/fc(To). Again, this quan- 
tity is extremely useful in kinetics because it can help to 
discern between plausible mechanisms of a reaction and 
because the ratio is easier to measure accurately than the 
rate constant itself. 

In the framework of the QI approximation, the T de- 
pendence can be evaluated as [22|] 



kqm ^ QrWo)AH{l3o) Cdd(/3) ^TlS) 

C'ddl.PoJ 
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The flux and energy variance are computed as for the 
KIEs, but the ratios of the partition functions and of the 
delta-delta correlation functions use a different type of 
TI, namely a TI with respect to the inverse temperature 
/3 [13, [131 . Taking advantage of the relations 



dlogQrW) 

d/3 

d log Cdd (/3) 
d/3 



and 
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To compute the EIE, ratios of partition functions cor- 
responding to different isotopologs of reactants and prod- 
ucts have to be computed. Unfortunately, higher level ab 
initio electronic structure methods are usually too ex- 
pensive to be used in the PI calculation directly. On the 
other hand, semiempirical or force field methods gener- 
ally do not achieve comparable accuracy. To take advan- 
tage of both the accuracy of ab initio methods and the 
rigor of PI treatment of nuclear motion, we combine the 
value of the EIE obtained in the harmonic approximation 
(HA) using a higher level method with the anharmonicity 
correction, rigorously computed with the PIMD method, 
using a lower level method. The anharmonicity correc- 
tion is calculated as 



(17) 



where AFpj'^ and AF^^ are the free energy differences 
computed with the PI method and the HA, respectively. 
The HA value of AF'^°'^ is obtained by Boltzmann aver- 
aging over all possible distinguishable conformations. 



AF{if = -ksTln 



«.E;!ri{exp(^)E;=iQ?:i5} 



(18) 

where Nj. is the number of "geometrically different iso- 
mers" of the reactant. By geometrically different iso- 
mers we mean species differing in their geometry, not 



4 



species differing only in positions of isotopically substi- 
tuted atoms. Ef is the electronic energy (including nu- 
clear repulsion) of the z*^ isomer, s,. is the symmetry 
factor, which can arise from the change of the wave func- 
tion due to the effects of indistinguishability of particles 
during isotopic substitution. Finally, Qr^ij are partition 
functions of the nuclear motion of Sr isotopomers. Np, 
Sp, Qp^ij denote analogous quantities for the product. 

Since in our approximation the electronic function 
does not change after isotopic substitution, the EIE 
is dominated by vibrational contributions. Therefore, 
we searched for an optimal electronic structure method 
among those studied by Merrick et al. [11], who tested 
the performance of several higher level methods by com- 
parison of the computed vibrational frequencies in the 
HA with experimental data for a set of 39 molecules. 
The B98/6-311-|-(2df,p) method, which we chose, has the 
root mean square error (RMSE) of the zero point energy 
(ZPE) equal to 0.31 kJ • mor^ and RMSEs of frequen- 
cies 31 cm^^. The lower level method used to compute 
the PI anharmonicity corrections was the AMI semiem- 
pirical method [29^, which, in the HA, reproduces the 
B98/6-311-h(2df,p) results very weh. 

Since the value of the EIE is expected to be close to 
unity, we have chosen a relatively high number of time 
slices P = 40 in the discretized path integral. The TI 
was performed with Simpson's rule using five values of 
A. As the dependence of the PIMD average of dF/dX is 
almost linear over the full range of A, five points were 
sufficient to get a converged result for the integral. For 
further details about the method, see Ref. [H. 



B. Kinetic isotope effects 

The ratio of partition functions needed for the KIE was 
computed in the same manner as the ratio for the EIE, 
with the exception that the number of imaginary time 
slices was set to P = 24. Other terms in the Eq. (fTTI) were 
computed with the PIMC method using the empirical 
valence bond (EVB) potential [13, [3l|, which allows the 
molecular mechanics potential to be fitted to match the 
B98/6-311-|-(2df,p) reaction barrier and Hessian at the 
transition state (33 |. 

The EVB potential is computed as 



^EVB = ^ {Vn + V22) - y \Vu\' \ j (19) 

where Vn and V22 are the molecular mechanics poten- 
tial energies of the reactant and product, respectively. 
These are the diabatic potentials and the diagonal terms 
in the symmetric 2x2 EVB matrix. In our calculations, 
they were obtained with the general AMBER force field 
(GAFF) [13. The off-diagonal term V12, i.e., the cou- 
pling between the two diabatic states, was calculated ac- 



cording to the formula of Schlegel and Sonnenberg [s^l : 

Fi2(r) = A[l + B-Ar+Ar-(C-|-aI)-Ar]exp(-a|Ar|V2). 

(20) 

The constants A, B, and C are chosen to match the 
barrier height and the Hessian of the ab initio potential 
at the transition state. 



C. Temperature dependence of the rate constant 

A similar procedure as for the KIE was used to com- 
pute the temperature dependence of the rate constant. 
Since the major contribution to the temperature de- 
pendence of the rate constant is due to the reaction 
barrier height, we computed the barrier height also by 
the single point coupled clusters CCSD(T)/cc-pVTZ and 
CCSD/aug-cc-pVTZ calculations on B98/6-311-f-(2df,p) 
geometries. 

In the PI calculation of the temperature dependence of 
the rate constant, several approaches were used to com- 
pute the ratio of reactant partition functions at differ- 
ent temperatures. In the first approach, denoted "QI 
GAFF(PIMC)", the ratio of the reactant partition func- 
tions was computed directly in the PIMC simulation with 
the GAFF force field. This could possibly cause a diffi- 
culty because in the EVB method, the Hessian of poten- 
tial energy surface at the transition state is fitted to the 
much more accurate B98/6-311-|-(2df,p) Hessian. At this 
point, the fitted surface has B98/6-311-|-(2df,p) proper- 
ties in the HA, with the anharmonicity given by a non- 
linear combination of the GAFF force fields describing 
the reactant and product. Therefore, a systematic error 
could arise due to the different properties of the B98 and 
GAFF surfaces in the HA. Fortunately, the harmonic vi- 
brational free energy computed using the GAFF force 
field is only about 0.5 kcal/mol lower than the B98/6- 
311-|-(2df,p) value. Consequently, the results obtained 
with the QI GAFF(PIMC) method are almost equal to 
the "QI B98(HA) GAFF(PIMG) anharm." results, ob- 
tained by augmenting the B98/6-311-t-(2df,p) free ener- 
gies in the HA with the anharmonicity PIMC correc- 
tion computed with the GAFF force field. The "QI 
B98(HA)" method uses only plain B98/6-31H-(2df,p) 
free energies in the HA. Finally, probably the most accu- 
rate "QI B98(HA) AMI (PIMC) anharm." method uses 
the B98/6-311+(2df,p) free energies in the HA with the 
anharmonicity PIMC correction computed with the AMI 
semiempirical potential. 



Used software 

All PIMD calculations were performed in Amber 10 
[ll|. The ab initio and B98 calculations as well as AMI 
semiempirical calculations in the HA were done in Gaus- 
sian 03 revision EOl [s^. All PIMC calculations were 
done using a PIMC code developed by one of us. 
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FIG. 1: The [1,5] hydrogen shift reaction in {3Z)- 
(5,5,5-^H3)penta-l,3-diene (l-5,5,5-d3) and in {3Z)-{1,1- 
^H2)penta-l,3-diene (1-1,1-^2)- If all contributions except 
for those due to symmetry factors were neglected, one would 
obtain approximate equilibrium constants Ki — 3 for the first 
reaction step and K2 = 2 for the second step (in both cases). 



IV. RESULTS 

In this section, the path integral formalism is applied to 
study the [1,5] sigmatropic hydrogen shift reaction in the 
(3Z)-penta-l,3-diene. Two isotopologs of (32')-penta- 
1,3-diene are considered: tri-deuterated (3Z)-(5,5,5- 
^H3)penta-l,3-diene (1-5,5,5-^3) and di-deuterated {3Z)- 
(l,l-2H2)penta-l,3-d iene (1-1,1-^2) (see Fig. [!}. Both 
isotopologs were used by Roth and Konig in their exper- 
imental study of the KIE in the [1,5] sigmatropic hydro- 
gen shift reaction (36)]. The value of the KIE obtained 
by Roth and Konig was too high to be explained classi- 
cally implying that quantum effects such as tunneling are 
important in this reaction. Here, we examine both the 
EIE and KIE for this reaction as well as the temperature 
dependence of the rate constant in the range studied in 
the experimental work [36i] . 
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FIG. 2: The KIE on the first [1,5] hydrogen shift reaction 
from Fig. [1] The displayed KIE is defined precisely in Eq. 
(|2ip of the text. The curve marked as exp*^ denotes the raw 
experimental data from Ref. [s^ whereas exp'' denotes the 
data computed from an Arrhenius fit in Ref. i36i . 

the equilibrium ratios beyond the symmetry effects. The 
value of the EIE for the first reaction step is greater than 
unity, meaning that the equilibrium of the first step re- 
action is shifted toward the product side for the lighter 
isotopolog l-l,l-d2 (where hydrogen is transferred) in 
comparison with the heavier 1-5,5,5-^3 (where the trans- 
ferred atom is deuterium). In the second step, the EIE 
is smaller than unity so the products are favored for the 
heavier isotopolog. Also, in contrast to the first step, hy- 
drogen is transferred in the heavier isotopolog and deu- 
terium in the lighter one. 



B. Kinetic isotope effects 



A. Equilibrium isotope effects 

The EIE for both steps of the [1,5] sigmatropic hydro- 
gen shift reaction is computed according to the Eq. ([7]) 
and listed in Table ID We have already computed the 
equilibrium ratios for both isotopologs in Ref. [l^. For 
the convenience of the reader. Table U also shows these 
ratios. 

The precision of the experimental ratios of Roth 
and Konig is only to a single significant digit, so it 
reflects only the ratios caused by the symmetry ef- 
fects and does not allow for a rigorous comparison 
with the computed EIE. Nevertheless, in Ref. [HI we 
were able to use the same methodology to compare 
the theoretical and experimental equilibrium ratios for 
a related compound (2,4,6,7,9-pentamethyl-5-methylene- 
ll,lla-dihydro-127J-naphthacene). As the experimental 
data of Doering et al. [STj were more accurate, a more 
rigorous test was possible, and led to the conclusion that 
the current methodology can reproduce quantitatively 



Computed values of the KIE on the first reaction step 
are shown in Figure [5] and Table |lll This KIE is defined 
as the ratio 

KIE^ Ml W2 ^ 1-5,5-..) 

A:(l-5,5,5-d3 l-l,l,5-d3) 

As can be seen from the table, the calculated and ex- 
perimental data agree very well; the largest relative dif- 
ference is around 5 %. Comparison with the transition 
state theory (TST) values, computed using the B98/6- 
311+(2df,p) barrier height and the partition functions in 
the HA [for reactants computed according to Eq. (IT5|) ]. 
clearly demonstrates the importance of quantum effects 
on the KIE. (Note that in TST the value of the KIE is 
actually independent of the barrier height). The third 
column of the table contains the KIE computed in Ref. 
[TtI , where the same PI method was used, but the barrier 
height and the Hessian of the transition state were fitted 
to MP2/6-31g(d) ab initio data. In addition, the GAFF 
force field 33] was used not only for the EVB potential 
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EIE 


EIE J 


= Ki{d2)/Ki 




AMI (PIMD) 




1.07 


0.96 


B98 (HA) +AM1 (PIMD) anharm. 




1.12 


0.95 


Equilibrium fraction 


1-5, 5, 5-^3 


l-l,l,5-d3 


l-l,5,5-d3 


AMI (PIMD) 


0.103 


0.296 


0.601 


B98 (HA) + AMI (PIMD) anharm. 


0.104 


0.293 


0.602 


Equilibrium fraction 


l-l,l-d2 


l-5,5-d2 


l-l,5-d2 


AMI (PIMD) 


0.099 


0.305 


0.597 


B98 (HA) +AM1 (PIMD) anharm. 


0.097 


0.307 


0.596 



TABLE I: Values of the EIE on the first (Ki) and second {K2) step of the reactions from Fig. [T] at 478.45 K. Results 
denoted by AMI (PIMD) were obtained with a PIMD calculation using the AMI semiempirical method. Results denoted by 
B98 (HA) + AMI (PIMD) anharm. were obtained using the B98 method in the harmonic approximation together with the 
anharmonicity correction computed with the AMI method. For reference, the bottom six rows show the equilibrium fractions 
of all isotopomers occurring in the reactions from Fig. [1] 
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TABLE II: The KIE on the first [1,5] hydrogen shift reaction 
from Fig. [T] The displayed KIE is defined precisely in Eq. 
(|21|l of the text. The column denoted by exp'' contains the 
raw experimental data from Ref. IIB whereas the column exp^ 
contains the data computed from an Arrhenius fit in Ref. ISq . 



but also for the reactant partition functions. Our results 
and results from Ref. T? differ by at most 10% of the KIE 
value. The main improvement in comparison with Ref. 
ITtI occurs at higher temperatures. However, the observed 
difference is not only due to the use of different methods, 
but partially also due to the way the EVB potential sur- 
face is constructed. As discussed in Ref. [13, the results 
are weakly dependent on the value of a in Eq. ([20)1 , and 
even though a flat plateau is obtained for a broad range 
of a values, varying a still can change the KIE by about 
10%. Generally, the optimal value of the parameter a can 
be different for different ab initio methods. Here we have 
intentionally chosen the same value of a = 0.9 as in Ref. 
[TtI . As the final results are still fairly close and both agree 
well with the experimental value, this demonstrates the 
relative robustness of the method once a proper value of 
a is chosen. In order to achieve an even higher accuracy 
and to remove a certain level of arbitrariness connected 
with the parameter a, one can use, e.g., the distributed 
gaussian method [s^. In this method, the EVB potential 
is fitted to match the Hessian of an ab initio potential at 
several points in the transition state region. 



C. Temperature dependence of the rate constant 

In contrast to the KIE which is virtually independent 
of the barrier height, the rate constant depends on the 
barrier height exponentially. Very accurate value of the 
barrier height is therefore essential for the calculation of 
this quantity or its temperature dependence. Table Hill 
contains barrier heights, relative energies of trans and 
gauche isomers of (3Z)-penta-l,3-diene, and TST val- 
ues of the rate constant kj^"^ for the 1-1,1-^2 ^ 1-5,5- 
d2, hydrogen shift computed at three different levels of 
theory using B98/6-311+(2df,p) optimized geometries. 
As can be seen from the table, the B98/6-311-t-(2df,p) 
barrier is the lowest one, in accordance with the well 
known fact that barrier heights are usually underesti- 
mated by the majority of currently used density func- 
tional methods (including B98) [38|. Besides B98, two 
other and generally more accurate methods used were the 
CCSD(T) method with the cc-pVTZ basis set and the 
CCSD method without triples correction but with the 
substantially larger aug-cc-pVTZ basis set. The differ- 
ence between these two methods is still almost 3 kcal/mol 
but at this point it is hard to decide which method is ac- 
tually closer to the real value of the barrier. 

The TST rate constants hsted in Table HIT] were com- 
puted using the B98/6-311-)-(2df,p) partition functions in 
the HA together with the barrier height of the method 
considered. Neglecting the effect of barrier recrossing 
and the anharmonicity of the potential energy surface, 
the TST rate constant should be smaller than the ex- 
perimental one, because the TST ignores the tunneling 
contribution that is quite important in this reaction (as 
demonstrated by the large value of the KIE) . As a conse- 
quence, comparison of kj^"^ with the experimental rate 
constant suggests that both the CCSD(T)/cc-pVTZ and 
especially the B98/6-311-|-(2df,p) methods give too low 
barriers. The TST rate constant computed with the 
CCSD/aug-cc-pVTZ barrier, on the other hand, agrees 
quite well with the theoretical expectations. 

The QI method provides several theoretical improve- 
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B98/6-311+(2df,p) CCSD(T)/cc 


pVTZ CCSD/aug-cc-pVTZ exp" 


Ae* (kcal/mol) 


37.3 39.9 


42.8 


Aef-*(kcal/mol) 


3.6 3.0 


2.0 


kj^'^ ■ 10^^ (s"') [478.45 K] 


364.6 27.8 


2.8 7.1 



"Ref. m 

TABLE III: The barrier height Ae* = — e*''''"'' measured from the trans conformer (3Z)-penta-l,3-diene which is the global 
minimum, the energy difference Ae^^* — e^^"'^'"^ — e''^^"^ between trans and gauche conformers, and the TST rate constant 
kj^'^ of the [1,5] hydrogen shift reaction l-l,l-(i2 — )■ 1-5,5-^2- 



1 

0.9 
0.8 
0.7 




\ Ql B98(HA) AM1(PIMC) anharm. i — i — i 
\ Ql B98(HA) GAFF(PIMC) anharm. 
\. Ql B98(HA) ' -1- - ■ 
\ QIGAFF(PiMC) ^ * - 
\» TST 


0.6 




■U, exp ■ 


0.5 






0.4 
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lOOO/nK) 

FIG. 3: Temperature dependence of the rate constant of the 
[1,5] hydrogen shift reaction l-l,l-(i2 — >■ l-5,5-d2 using the 
CCSD/aug-cc-pVTZ barrier height. The y axis is plotted in 
the logarithmic scale. The description of methods used can 
be found in Computational details. 

ments over the TST: Most importantly, unlike the TST, 
the Ql approximation includes tunneling and other quan- 
tum effects. Moreover, its PI implementation facili- 
tates the calculation of the anharmonicity effects on the 
partition functions and delta-delta correlation functions. 
Methods used to include the anharmonicity correction to 
the B98/6-311-|-(2df,p) harmonic ratios are described in 
the Computational details. The temperature dependence 
of the Ql rate constant with the CCSD/aug-cc-pVTZ 
barrier is shown in Figure [31 

As can be seen in Figure [Sj both the Ql and TST re- 
sults agree well with the relatively broad and imprecise 
experimental temperature dependence. Neglecting again 
the recrossing and anharmonicity effects, one can expect 
that the value of kj^'^ will be smaller than the value of 
k'^^, mainly due to the tunneling contribution to the rate 

that is captured by k'^^ but ignored in kj^'^ . Moreover, 
as the relative importance of the tunneling increases with 
the decreasing temperature, the temperature dependence 
of the k'^^ should be weaker than the dependence of kj^^ . 
Figure [3] shows that this is indeed true. (The tempera- 
ture dependence of the relative importance of tunneling 
and of other nuclear quantum effects can be sometimes 
more complicated due to other effects such as the tem- 



perature dependence of the distance between the donor 
and acceptor [2lj.) 

Another interesting feature visible in the graph is the 
importance of the anharmonicity correction, which shows 
that the change in the temperature dependence caused 
by the anharmonicity can be equivalent to the change of 
the reaction barrier by several kcal/mol. It is reassuring 
that the more accurate AMI correction influences the re- 
sult more than the correction computed using the GAFF 
force field. This observation agrees well with the expec- 
tations based on the fact that in the GAFF force field 
the bond stretching coordinates are actually harmonic 
and the anharmonicity is caused only by the other terms 
of the force field. Nevertheless, both corrections decrease 
the HA free energy of nuclear motion through the low- 
ering of the ZPE. Finally, the curves of the relative rate 
constant in the Arrhenius plot are almost exactly straight 
lines, irrespective of the theoretical method used. This 
suggests that it would be very hard to observe any non- 
Arrhenius behavior in this small temperature range even 
with a far more precise experimental setup. 

The temperature dependence of fc^j computed with the 
CCSD(T)/cc-pVTZ barrier height is shown in Figure HI 
With this barrier, the temperature dependence of the 
more accurate Ql method is only very slightly weaker 
than the experimental dependence, whereas the TST rate 
dependence is still within the range of the experimental 
values. However, this apparently excellent agreement of 
the conventional TST with the experiment is most likely 
due to the cancellation of the tunneling and anharmonic- 
ity corrections. All the other characteristics of the graph 
are the same as with the CCSD/aug-cc-pVTZ barrier. 
When the lowest and probably the least accurate B98 / 6- 
311-|-(2df,p) barrier is used, curves for all methods stay 
above the experimental data in the Arrhenius plot. 

V. CONCLUSIONS 

We have applied a general path integral methodology 
for computing the EIE, KIE, and the temperature depen- 
dence of the rate constant to the [1,5] sigmatropic hydro- 
gen shift in pentadiene. In case of the EIE, the computed 
result has a higher precision than the experiment [36| . 
therefore our result can be considered as a prediction of 
the deviation of the exact EIE from the symmetry deter- 
mined EIE. In case of the KIE, the accuracy of the previ- 
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lOOO/nK) 

FIG. 4: Temperature dependence of the rate constant of the 
[1,5] hydrogen shift reaction 1-1,1-^2 — >■ l-5,5-(i2 using the 
CCSD(T)/cc-pVTZ barrier height. The y axis is plotted in 
the logarithmic scale. The description of methods used can 
be found in Computational details. 

ous result 17] was improved by using a combination of a 
high level electronic structure calculation within the har- 
monic approximation with a path integral anharmonicity 
correction using a lower level method. Finally, the result 
for the temperature dependence of the rate constant is 
the first application of the methodology from Ref. 22 
to a molecule with more than three atoms. As for the 



KIE, the temperature dependence confirms the impor- 
tance of tunneling and anharmonicity effects. However, 
unlike for the KIE, the accuracy of the energy barrier 
plays an important role and therefore a high level method 
is required. According to our results, the CCSD/aug-cc- 
pVDZ barrier height of 42.8 kcal/mol seems to be the 
most accurate. 

While we have considered a gas phase molecule with 13 
atoms, significantly larger systems can be treated by the 
presented PI approach since only the light atoms par- 
ticipating in the reaction require a large number P of 
replicas. Heavier atoms or atoms less important in a 
given problem (e.g., atoms far from the active site of 
an enzyme) can be treated either with a small value of 
P or fully classically (i.e., with P — 1). The efficiency 
will then be dominated by the cost of an accurate po- 
tential energy for which various strategies are known al- 
ready from classical molecular dynamics. Among these 
are the recent implementations (see, e.g., Ref. of po- 
larizable force fields (i^, |4l| or the accelerated sampling 
using hybrid quantum-mechanical / molecular-mechanical 
potentials [42]. A combination of such techniques with 
the methodology described in this paper will be the sub- 
ject of future applications. 
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